%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  user function
%  1
%
function [f1,f2] = rhs(x,y,Re, caseNum)
    nu = 1/Re;
    
    switch caseNum
        case 0  % stream function with analytic expression
    x2 = x.*x;    y2 = y.*y;
    x3 = x.*x2;   y3 = y.*y2;
    x4 = x2.*x2;  y4 = y2.*y2;
    x5 = x2.*x3;  y5 = y2.*y3;
    x6 = x3.*x3;  y6 = y3.*y3;
    x7 = x3.*x4;  y7 = y3.*y4;
    f1 = 4*x3.*y2 - 16*x3.*y3 - 20*x4.*y2 + 28*x3.*y4 + ...
        80*x4.*y3 + 36*x5.*y2 - 24*x3.*y5 - 140*x4.*y4 - ... 
        144*x5.*y3 - 28*x6.*y2 + 8*x3.*y6 + 120*x4.*y5 + ...
        252*x5.*y4 + 112*x6.*y3 + 8*x7.*y2 - 40*x4.*y6 - ...
        216*x5.*y5 - 196*x6.*y4 - 32*x7.*y3 + 72*x5.*y6 + ...
        168*x6.*y5 + 56*x7.*y4 - 56*x6.*y6 - 48*x7.*y5 + ...
        16*x7.*y6 - 4*nu*y + 12*nu*x2 - 24*nu*x3 + 12*nu*x4 + ...
        12*nu*y2 - 8*nu*y3 + 3*x2 - 72*nu*x.*y2 - 48*nu*x2.*y + ...
        48*nu*x.*y3 + 48*nu*x3.*y - 24*nu*x4.*y + 72*nu*x2.*y2 - ...
        48*nu*x2.*y3 + 24*nu*x.*y;
    f2 = 4*x2.*y3 - 20*x2.*y4 - 16*x3.*y3 + 36*x2.*y5 + ...
        80*x3.*y4 + 28*x4.*y3 - 28*x2.*y6 - 144*x3.*y5 - ...
        140*x4.*y4 - 24*x5.*y3 + 8*x2.*y7 + 112*x3.*y6 + ...
        252*x4.*y5 + 120*x5.*y4 + 8*x6.*y3 - 32*x3.*y7 - ...
        196*x4.*y6 - 216*x5.*y5 - 40*x6.*y4 + 56*x4.*y7 + ...
        168*x5.*y6 + 72*x6.*y5 - 48*x5.*y7 - 56*x6.*y6 + ...
        16*x6.*y7 + 4*nu*x - 12*nu*x2 + 8*nu*x3 - 12*nu*y2 + ...
        24*nu*y3 - 12*nu*y4 + 3*y2 + 48*nu*x.*y2 + 72*nu*x2.*y - ...
        48*nu*x.*y3 - 48*nu*x3.*y + 24*nu*x.*y4 - 72*nu*x2.*y2 + ...
        48*nu*x3.*y2 - 24*nu*x.*y;
    
        case 1  % kovasznay's exact solution
            
    
        case 2  % square cavity
            f1 = 0*x;
            f2 = f1;
    end
            
end